C
      SUBROUTINE FMI3AL(INUHF,IOUT,TRNOP,NPERFL,ISS,IVER)
C **********************************************************************
C THIS SUBROUTINE CHECKS FLOW-TRANSPORT LINK FILE AND ALLOCATES SPACE
C FOR ARRAYS THAT MAY BE NEEDED BY FLOW MODEL-INTERFACE (FMI) PACKAGE.
C **********************************************************************
C LAST MODIFIED: 06-23-98
C

      IMPLICIT	NONE
      INTEGER	INUHF,IOUT,NPERFL,
     &		MTWEL,MTDRN,MTRCH,MTEVT,MTRIV,MTGHB,MTCHD,ISS,IVER
      CHARACTER VERSION*11
      LOGICAL	TRNOP(10),FWEL,FDRN,FRCH,FEVT,FRIV,FGHB
      COMMON /FC/FWEL,FDRN,FRCH,FEVT,FRIV,FGHB
	INTEGER*4 RESULT
C
C--PRINT PACKAGE NAME AND VERSION NUMBER
      WRITE(IOUT,1030) INUHF
 1030 FORMAT(1X,'FMI3 -- FLOW MODEL INTERFACE PACKAGE,',
     & ' VER DoD_3.0, JUNE 1998, INPUT READ FROM UNIT',I3)
C
C--INITIALIZE
      FWEL=.FALSE.
      FDRN=.FALSE.
      FRCH=.FALSE.
      FEVT=.FALSE.
      FRIV=.FALSE.
      FGHB=.FALSE.
      MTCHD=0
      ISS=1
      NPERFL=0
      IVER=1
C
C--READ HEADER OF FLOW-TRANSPORT LINK FILE
      READ(INUHF,ERR=500,END=500) VERSION,MTWEL,MTDRN,MTRCH,MTEVT,
     & MTRIV,MTGHB,MTCHD,ISS,NPERFL
C
C--DETERMINE WHICH FLOW COMPONENTS USED IN FLOW MODEL
      IF(VERSION(1:4).EQ.'MT3D') THEN
	IVER=2
	IF(MTWEL.GT.0) FWEL=.TRUE.
	IF(MTDRN.GT.0) FDRN=.TRUE.
	IF(MTRCH.GT.0) FRCH=.TRUE.
	IF(MTEVT.GT.0) FEVT=.TRUE.
	IF(MTRIV.GT.0) FRIV=.TRUE.
	IF(MTGHB.GT.0) FGHB=.TRUE.
      ENDIF
C
C--DETERMINE IF MT3DMS SSM PACKAGE IS REQUIRED
  200 IF(.NOT.TRNOP(3)) THEN
	IF(FWEL.OR.FDRN.OR.FRCH.OR.FEVT.OR.FRIV.OR.FGHB.OR.
     &	 MTCHD.GT.0.OR.ISS.EQ.0) THEN
	  WRITE(*,300)
C-------EMRL JIG
        call stopfile
C-------EMRL JIG
	  STOP
	ENDIF
      ENDIF
  300 FORMAT(/1X,'ERROR: THE RT3D [SSM] PACKAGE MUST BE USED',
     & ' FOR THE CURRENT SIMULATION',
     & /1X,'BECAUSE ONE OF THE FOLLOWING APPLIES:',
     & /1X,'1. THE FLOW MODEL INCLUDES A SINK/SOURCE PACKAGE;',
     & /1X,'2. THE FLOW MODEL CONTAINS CONSTANT-HEAD CELLS; OR',
     & /1X,'3. THE FLOW MODEL IS TRANSIENT.')
C
C--PRINT KEY INFORMATION OF THE FLOW MODEL
      IF(ISS.EQ.0) THEN
        WRITE(IOUT,310)
      ELSE
        WRITE(IOUT,320)
      ENDIF
      WRITE(IOUT,330)
      IF(FWEL) WRITE(IOUT,340)
      IF(FDRN) WRITE(IOUT,342)
      IF(FRCH) WRITE(IOUT,344)
      IF(FEVT) WRITE(IOUT,346)
      IF(FRIV) WRITE(IOUT,348)
      IF(FGHB) WRITE(IOUT,350)
      IF(MTCHD.GT.0) WRITE(IOUT,360)
      WRITE(IOUT,'(1X)')
  310 FORMAT(1X,'FLOW SOLUTION IS TRANSIENT')
  320 FORMAT(1X,'FLOW SOLUTION IS STEADY-STATE')
  330 FORMAT(1X,'SINK/SOURCE COMPONENTS INCLUDED IN FLOW MODEL:')
  340 FORMAT(1X,'o WELL')
  342 FORMAT(1X,'o DRAIN')
  344 FORMAT(1X,'o RECHARGE')
  346 FORMAT(1X,'o EVAPOTRANSPIRATION')
  348 FORMAT(1X,'o RIVER OR STREAM')
  350 FORMAT(1X,'o GENERAL-HEAD-DEPENDENT BOUNDARY')
  360 FORMAT(1X,'o CONSTANT-HEAD CELLS')
C
C--DONE, RETURN
      GOTO 1000
C
C--ISSUE A WARNING IF THE LINK FILE IS SAVED BY VER 1 OF LKMT PACKAGE
  500 REWIND(INUHF)
      WRITE(IOUT,505)
      WRITE(*,510)
  505 FORMAT(1X,'FLOW-TRANSPORT LINK FILE SAVED BY',
     & ' VER 1 OF [LKMT] PACKAGE'/)
  510 FORMAT(/1X,'WARNING: FLOW-TRANSPORT LINK FILE SAVED BY',
     & ' VER 1 OF [LKMT] PACKAGE;'/1X,9X,'RUN FLOW MODEL WITH',
     & ' VER 2 OR ABOVE OF [LKMT] PACKAGE',
     & /1X,9X,'TO TAKE FULL ADVANTAGE OF NEW CAPABILITIES IN RT3D.')
C
 1000 RETURN
      END
C
C
      SUBROUTINE FMI3RP1(INUF,IOUT,KPER,KSTP,NCOL,NROW,NLAY,NCOMP,
     & FPRT,LAYCON,ICBUND,HORIGN,DH,PRSITY,DELR,DELC,DZ,XBC,YBC,ZBC,
     & QSTO,COLD,CNEW,RETA,QX,QY,QZ,DTRACK,DTRACK2,THKMIN,ISS,IVER)
C **********************************************************************
C THIS SUBROUTINE READS SATURATED CELL THICKNESS, FLUXES ACROSS CELL
C INTERFACES, AND FLOW RATE TO OR FROM TRANSIENT STORAGE
C FROM AN UNFORMATTED FILE SAVED BY THE FLOW MODEL, AND PREPARES THEM
C IN THE FORMS NEEDED BY THE TRANSPORT MODEL.
C **********************************************************************
C last modified: 06-23-98
C

      IMPLICIT	NONE
      INTEGER	INUF,IOUT,NCOL,NROW,NLAY,LAYCON,ICBUND,J,I,K,KPER,KSTP,
     &          JTRACK,ITRACK,KTRACK,IVER,ISS,NCOMP,INDEX
      REAL	DH,PRSITY,DELR,DELC,DZ,XBC,YBC,ZBC,HORIGN,QX,QY,QZ,WW,
     &		WTBL,THKSAT,DTRACK,TK,CNEW,COLD,RETA,CTMP,CREWET,QSTO,
     &		THKMIN,THKMIN0,DTRACK2
      CHARACTER FPRT*1,TEXT*16
      DIMENSION LAYCON(NLAY),ICBUND(NCOL,NROW,NLAY,NCOMP),
     &		DH(NCOL,NROW,NLAY),
     &		DELR(NCOL),DELC(NROW),DZ(NCOL,NROW,NLAY),XBC(NCOL),
     &		YBC(NROW),ZBC(NCOL,NROW,NLAY),QX(NCOL,NROW,NLAY),
     &		QY(NCOL,NROW,NLAY),QZ(NCOL,NROW,NLAY),
     &		PRSITY(NCOL,NROW,NLAY),CNEW(NCOL,NROW,NLAY,NCOMP),
     &		COLD(NCOL,NROW,NLAY,NCOMP),RETA(NCOL,NROW,NLAY,NCOMP),
     &		QSTO(NCOL,NROW,NLAY)
	INTEGER*4 RESULT
C
C--READ SATURATED THICKNESS (UNIT: L).
      IF(IVER.EQ.2) THEN
	TEXT='THKSAT'
      ELSEIF(IVER.EQ.1) THEN
	TEXT='HEAD'
      ENDIF
      CALL READHQ(INUF,IOUT,NCOL,NROW,NLAY,KSTP,KPER,TEXT,DH,FPRT)
C
C--READ RIGHT-FACE FLOW TERMS IF MORE THAN ONE COLUMN (UNIT: L**3/T).
      IF(NCOL.LT.2) GOTO 100
      TEXT='QXX'
      CALL READHQ(INUF,IOUT,NCOL,NROW,NLAY,KSTP,KPER,TEXT,QX,FPRT)
C
C--READ FRONT-FACE FLOW TERMS IF MORE THAN ONE ROW (UNIT: L**3/T).
  100 IF(NROW.LT.2) GOTO 110
      TEXT='QYY'
      CALL READHQ(INUF,IOUT,NCOL,NROW,NLAY,KSTP,KPER,TEXT,QY,FPRT)
C
C--READ LOWER-FACE FLOW TERMS IF MORE THAN ONE LAYER (UNIT: L**3/T).
  110 IF(NLAY.LT.2) GOTO 120
      TEXT='QZZ'
      CALL READHQ(INUF,IOUT,NCOL,NROW,NLAY,KSTP,KPER,TEXT,QZ,FPRT)
C
C--READ STORAGE TERM (UNIT: L**3/T).
  120 TEXT='STO'
      IF(IVER.EQ.2.AND.ISS.EQ.0) THEN
	CALL READHQ(INUF,IOUT,NCOL,NROW,NLAY,KSTP,KPER,TEXT,QSTO,FPRT)
      ENDIF
C
C--SET ICBUND=0 IF CELL IS DRY OR INACTIVE (H=1.E30)
C--AND REACTIVATE DRY CELL IF REWET AND ASSIGN CONC AT REWET CELL
      DO K=1,NLAY
	DO I=1,NROW
	  DO J=1,NCOL
	    IF(ABS(DH(J,I,K)-1.E30).LT.1.E-5) THEN
	      ICBUND(J,I,K,1)=0
	    ELSEIF(ICBUND(J,I,K,1).EQ.0.AND.PRSITY(J,I,K).GT.0) THEN
	      ICBUND(J,I,K,1)=30000
	      DO INDEX=1,NCOMP
		CTMP=CREWET(NCOL,NROW,NLAY,CNEW(1,1,1,INDEX),
     &			    ICBUND,XBC,YBC,ZBC,J,I,K)
		CTMP=(COLD(J,I,K,INDEX)*(RETA(J,I,K,INDEX)-1.0)+CTMP)
     &		    /RETA(J,I,K,INDEX)
		CNEW(J,I,K,INDEX)=CTMP
	      ENDDO
	    ENDIF
	  ENDDO
	ENDDO
      ENDDO
C
C--GET SATURATED THICKNESS DZ FOR CONFINED LAYERS
C--FOR FILE SAVED BY LKMT PACKAGE VERSION 2
      IF(IVER.EQ.2) THEN
C
	DO K=1,NLAY
	  DO I=1,NROW
	    DO J=1,NCOL
	      IF(ICBUND(J,I,K,1).NE.0) THEN
		IF(LAYCON(K).EQ.0 .OR. INT(DH(J,I,K)).EQ.-111) THEN
		  DH(J,I,K)=DZ(J,I,K)
		ENDIF
                THKMIN0=THKMIN*DZ(J,I,K)
                IF((DH(J,I,K).LT.THKMIN0).or.(DH(J,I,K).eq.0.0)) THEN
                  WRITE(IOUT,105) DH(J,I,K),K,I,J,THKMIN0
		  ICBUND(J,I,K,1)=0
		ENDIF
	      ENDIF
	    ENDDO
	  ENDDO
	ENDDO
C
C--FOR FILE SAVED BY LKMT PACKAGE VERSION 1
      ELSEIF(IVER.EQ.1) THEN
C
	DO K=1,NLAY
	  DO I=1,NROW
	    DO J=1,NCOL
C
	      IF(ICBUND(J,I,K,1).NE.0) THEN
		IF(LAYCON(K).EQ.0) THEN
		  THKSAT=DZ(J,I,K)
		ELSE
		  WTBL=HORIGN-DH(J,I,K)
		  THKSAT=ZBC(J,I,K)+0.5*DZ(J,I,K)-WTBL
		  THKSAT=MIN(THKSAT,DZ(J,I,K))
		ENDIF
		IF(THKSAT.LE.0) THEN
		  WRITE(*,355) K,I,J
C-------EMRL JIG
            call stopfile
C-------EMRL JIG
		  STOP
		ELSE
		  DH(J,I,K)=THKSAT
		ENDIF
                THKMIN0=THKMIN*DZ(J,I,K)
                IF((DH(J,I,K).LT.THKMIN0).or.(DH(J,I,K).eq.0.0)) THEN
                  WRITE(IOUT,105) DH(J,I,K),K,I,J,THKMIN0
		  ICBUND(J,I,K,1)=0
		ENDIF
	      ENDIF
C
	    ENDDO
	  ENDDO
	ENDDO
C
      ENDIF
C
  105 FORMAT(/1X,'SATURATED THICKNESS =',G13.5,'AT CELL (K,I,J):',3I5,
     & /1X,'BELOW SPECIFIED MINIMUM =',G13.5,'CELL RESET AS INACTIVE')
  355 FORMAT(/1X,'ERROR: SATURATED THICKNESS CALCULATED TO BE ZERO',
     &	' OR NEGATIVE FOR ACTIVE CELL'/1X,'       AT LAYER=',I2,
     &	'  ROW=',I4,'  COLUMN=',I4,
     &	/1X,'CHECK [HTOP] AND [DZ] ARRAYS IN [BTN] INPUT FILE.')
C
C--DETERMINE MAXIMUM TIME INCREMENT DURING WHICH ANY PARTICLE
C--CANNOT MOVE MORE THAN ONE CELL IN ANY DIRECTION.
      DTRACK=1.E30
C
      IF(NCOL.LT.2) GOTO 410
      DO K=1,NLAY
	DO I=1,NROW
	  DO J=2,NCOL
C
	    IF(ICBUND(J,I,K,1).NE.0) THEN
	      TK=0.5*(QX(J-1,I,K)+QX(J,I,K))
	      IF(TK.EQ.0) CYCLE
	      TK=DELR(J)*DELC(I)*DH(J,I,K)*PRSITY(J,I,K)/TK
	      IF(ABS(TK).LT.DTRACK) THEN
		DTRACK=ABS(TK)
		JTRACK=J
		ITRACK=I
		KTRACK=K
	      ENDIF
	    ENDIF
C
	  ENDDO
	ENDDO
      ENDDO
C
  410 IF(NROW.LT.2) GOTO 420
      DO K=1,NLAY
	DO J=1,NCOL
	  DO I=2,NROW
C
	    IF(ICBUND(J,I,K,1).NE.0) THEN
	      TK=0.5*(QY(J,I-1,K)+QY(J,I,K))
	      IF(TK.EQ.0) CYCLE
	      TK=DELR(J)*DELC(I)*DH(J,I,K)*PRSITY(J,I,K)/TK
	      IF(ABS(TK).LT.DTRACK) THEN
		DTRACK=ABS(TK)
		JTRACK=J
		ITRACK=I
		KTRACK=K
	      ENDIF
	    ENDIF
C
	  ENDDO
	ENDDO
      ENDDO
C
  420 IF(NLAY.LT.2) GOTO 430
      DO J=1,NCOL
	DO I=1,NROW
	  DO K=2,NLAY
C
	    IF(ICBUND(J,I,K,1).NE.0) THEN
	      TK=0.5*(QZ(J,I,K-1)+QZ(J,I,K))
	      IF(TK.EQ.0) CYCLE
	      TK=DELR(J)*DELC(I)*DH(J,I,K)*PRSITY(J,I,K)/TK
	      IF(ABS(TK).LT.DTRACK) THEN
		DTRACK=ABS(TK)
		JTRACK=J
		ITRACK=I
		KTRACK=K
	      ENDIF
	    ENDIF
C
	  ENDDO
	ENDDO
      ENDDO
C
C--PRINT INFORMATION ON DTRACK
  430 WRITE(IOUT,500) DTRACK,KTRACK,ITRACK,JTRACK
  500 FORMAT(/1X,'MAXIMUM STEPSIZE DURING WHICH ANY PARTICLE CANNOT',
     & ' MOVE MORE THAN ONE CELL'/1X,'=',G11.4,
     & '(WHEN MIN. R.F.=1)  AT K=',I4,', I=',I4,
     & ', J=',I4)
C
C--DETERMINE STABILITY CRITERION ASSOCIATED WITH EXPLICIT FINITE
C--DIFFERENCE SOLUTION OF THE ADVECTION TERM
      DTRACK2=1.E30
C
      DO K=1,NLAY
	DO I=1,NROW
	  DO J=1,NCOL
	    IF(ICBUND(J,I,K,1).EQ.0) CYCLE
	    TK=0.
	    IF(J.GT.1) TK=TK+MAX(ABS(QX(J-1,I,K)),ABS(QX(J,I,K)))
	    IF(I.GT.1) TK=TK+MAX(ABS(QY(J,I-1,K)),ABS(QY(J,I,K)))
	    IF(K.GT.1) TK=TK+MAX(ABS(QZ(J,I,K-1)),ABS(QZ(J,I,K)))
	    IF(TK.EQ.0) CYCLE
	    TK=DELR(J)*DELC(I)*DH(J,I,K)*PRSITY(J,I,K)/TK
	    IF(TK.LT.DTRACK2) THEN
	      DTRACK2=TK
	      JTRACK=J
	      ITRACK=I
	      KTRACK=K
	    ENDIF
	  ENDDO
	ENDDO
      ENDDO
C
C--PRINT INFORMATION ON DTRACK2
      WRITE(IOUT,550) DTRACK2,KTRACK,ITRACK,JTRACK
  550 FORMAT(/1X,'MAXIMUM STEPSIZE WHICH MEETS STABILITY CRITERION',
     & ' OF THE ADVECTION TERM'/1X,
     & '(FOR PURE FINITE-DIFFERENCE OPTION, MIXELM=0) '/1X,'=',G11.4,
     & '(WHEN MIN. R.F.=1)  AT K=',I4,', I=',I4,', J=',I4)
C
C--DIVIDE VOLUMETRIC QX, QY AND QZ BY AREAS
C--TO GET SPECIFIC DISCHAGES ACROSS EACH CELL INTERFACE
      IF(NCOL.LT.2) GOTO 910
      DO K=1,NLAY
	DO I=1,NROW
	  DO J=1,NCOL-1
	    WW=DELR(J+1)/(DELR(J+1)+DELR(J))
	    THKSAT=DH(J,I,K)*WW+DH(J+1,I,K)*(1.-WW)
	    IF(THKSAT.LE.0.OR.ICBUND(J,I,K,1).EQ.0) THEN
	      QX(J,I,K)=0
              IF(J.GT.1) QX(J-1,I,K)=0.
            ELSE
	      QX(J,I,K)=QX(J,I,K)/(DELC(I)*THKSAT)
	    ENDIF
	  ENDDO
          IF(ICBUND(NCOL,I,K,1).EQ.0) QX(NCOL-1,I,K)=0.
        ENDDO
      ENDDO
C
  910 IF(NROW.LT.2) GOTO 920
      DO K=1,NLAY
	DO J=1,NCOL
	  DO I=1,NROW-1
	    WW=DELC(I+1)/(DELC(I+1)+DELC(I))
	    THKSAT=DH(J,I,K)*WW+DH(J,I+1,K)*(1.-WW)
	    IF(THKSAT.LE.0.OR.ICBUND(J,I,K,1).EQ.0) THEN
	      QY(J,I,K)=0
              IF(I.GT.1) QY(J,I-1,K)=0.
            ELSE
	      QY(J,I,K)=QY(J,I,K)/(DELR(J)*THKSAT)
	    ENDIF
	  ENDDO
          IF(ICBUND(J,NROW,K,1).EQ.0) QY(J,NROW-1,K)=0.
        ENDDO
      ENDDO
C
  920 IF(NLAY.LT.2) GOTO 990
      DO J=1,NCOL
	DO I=1,NROW
          DO K=1,NLAY
	    THKSAT=DH(J,I,K)
	    IF(THKSAT.LE.0.OR.ICBUND(J,I,K,1).EQ.0) THEN
	      QZ(J,I,K)=0
              IF(K.GT.1) QZ(J,I,K-1)=0.
            ELSE
	      QZ(J,I,K)=QZ(J,I,K)/(DELR(J)*DELC(I))
	    ENDIF
	  ENDDO
	ENDDO
      ENDDO
C
C--DIVIDE STORAGE BY CELL VOLUME TO GET DIMENSION (1/TIME)
  990 CONTINUE
      DO K=1,NLAY
	DO I=1,NROW
	  DO J=1,NCOL
	    THKSAT=DH(J,I,K)
	    IF(THKSAT.LE.0.OR.ICBUND(J,I,K,1).EQ.0) THEN
	      QSTO(J,I,K)=0
	    ELSE
	      QSTO(J,I,K)=QSTO(J,I,K)/(THKSAT*DELR(J)*DELC(I))
	    ENDIF
	  ENDDO
	ENDDO
      ENDDO
C
C--SYNCHRONIZE ICBUND CONDITIONS OF ALL SPECIES
      IF(NCOMP.EQ.1) GOTO 999
      DO K=1,NLAY
	DO I=1,NROW
	  DO J=1,NCOL
	    DO INDEX=2,NCOMP
	      IF(ICBUND(J,I,K,INDEX).GE.0) THEN
                ICBUND(J,I,K,INDEX)=IABS(ICBUND(J,I,K,1))
	      ELSEIF(ICBUND(J,I,K,1).EQ.0) THEN
		ICBUND(J,I,K,INDEX)=0
	      ENDIF
	    ENDDO
	  ENDDO
	ENDDO
      ENDDO
C
C--RETURN
  999 RETURN
      END
C
C
      SUBROUTINE FMI3RP2(INUF,IOUT,KPER,KSTP,NCOL,NROW,NLAY,NCOMP,
     & FPRT,LAYCON,ICBUND,DH,PRSITY,DELR,DELC,
     & IRCH,RECH,IEVT,EVTR,MXSS,NSS,NTSS,SS,BUFF,DTSSM)
C **********************************************************************
C THIS SUBROUTINE READS THE LOCATIONS AND FLOW RATES OF SINKS & SOURCES
C FROM AN UNFORMATTED FILE SAVED BY THE FLOW MODEL, AND PREPARES THEM
C IN THE FORMS NEEDED BY THE TRANSPORT MODEL.
C **********************************************************************
C last modified: 06-23-98
C

      IMPLICIT	NONE
      INTEGER	INUF,IOUT,NCOL,NROW,NLAY,MXSS,NSS,LAYCON,ICBUND,J,I,K,
     &		NUM,KPER,KSTP,IRCH,IEVT,NTSS,IQ,KSSM,ISSM,JSSM,
     &          JJ,II,KK,JM1,JP1,IM1,IP1,KM1,KP1,NCOMP,INDEX
      REAL	DH,PRSITY,DELR,DELC,SS,BUFF,VOLAQU,EVTR,RECH,DTSSM,TM
      LOGICAL	FWEL,FDRN,FRCH,FEVT,FRIV,FGHB
      CHARACTER FPRT*1,TEXT*16
      DIMENSION LAYCON(NLAY),ICBUND(NCOL,NROW,NLAY,NCOMP),
     &		DH(NCOL,NROW,NLAY),
     &		DELR(NCOL),DELC(NROW),PRSITY(NCOL,NROW,NLAY),
     &		IRCH(NCOL,NROW),RECH(NCOL,NROW),IEVT(NCOL,NROW),
     &		EVTR(NCOL,NROW),SS(6,MXSS),BUFF(NCOL,NROW,NLAY)
      COMMON /FC/FWEL,FDRN,FRCH,FEVT,FRIV,FGHB
	INTEGER*4 RESULT
C
C--RESET NUMBER OF POINT SOURCES OF SPECIFIED CONCENTRATIONS [NSS]
C--AND MAKE [ICBUND] VALUE AT ACTIVE CELLS EQUAL TO 1
      NTSS=NSS
      DO K=1,NLAY
	DO I=1,NROW
	  DO J=1,NCOL
	    IF(ICBUND(J,I,K,1).GT.0) ICBUND(J,I,K,1)=1
	  ENDDO
	ENDDO
      ENDDO
C
C--READ CONSTANT-HEAD FLOW TERM (UNIT: L**3/T).
      TEXT='CNH'
      IQ=1
      CALL READPS(INUF,IOUT,NCOL,NROW,NLAY,KSTP,KPER,TEXT,
     & BUFF,IQ,MXSS,NTSS,NSS,SS,ICBUND,FPRT)
C
C--READ WELL FLOW TERM (L**3/T) IF WELL OPTION USED IN FLOW MODEL.
  230 IF(.NOT.FWEL) GOTO 300
      TEXT='WEL'
      IQ=2
      CALL READPS(INUF,IOUT,NCOL,NROW,NLAY,KSTP,KPER,TEXT,
     & BUFF,IQ,MXSS,NTSS,NSS,SS,ICBUND,FPRT)
C
C--READ DRAIN FLOW TERM (L**3/T) IF DRAIN OPTION USED IN FLOW MODEL.
  300 IF(.NOT.FDRN) GOTO 400
      TEXT='DRN'
      IQ=3
      CALL READPS(INUF,IOUT,NCOL,NROW,NLAY,KSTP,KPER,TEXT,
     & BUFF,IQ,MXSS,NTSS,NSS,SS,ICBUND,FPRT)
C
C--READ RECHARGE FLOW TERM (L**3/T)
C--IF RECHARGE OPTION USED IN FLOW MODEL
  400 IF(.NOT.FRCH) GOTO 500
      TEXT='RCH'
      CALL READDS(INUF,IOUT,NCOL,NROW,NLAY,KSTP,KPER,TEXT,
     & RECH,IRCH,FPRT)
C
C--READ E.T. FLOW TERM (L**3/T) IF E.T. OPTION USED IN FLOW MODEL
  500 IF(.NOT.FEVT) GOTO 600
      TEXT='EVT'
      CALL READDS(INUF,IOUT,NCOL,NROW,NLAY,KSTP,KPER,TEXT,
     & EVTR,IEVT,FPRT)
C
C--READ RIVER FLOW TERM (L**3/T) IF RIVER OPTION USED IN FLOW MODEL.
  600 IF(.NOT.FRIV) GOTO 700
      TEXT='RIV'
      IQ=4
      CALL READPS(INUF,IOUT,NCOL,NROW,NLAY,KSTP,KPER,TEXT,
     & BUFF,IQ,MXSS,NTSS,NSS,SS,ICBUND,FPRT)
C
C--READ GERENAL HEAD DEPENDENT BOUNDARY FLOW TERM (L**3/T)
C--IF GHB OPTION IS USED IN FLOW MODEL.
  700 IF(.NOT.FGHB) GOTO 800
      TEXT='GHB'
      IQ=5
      CALL READPS(INUF,IOUT,NCOL,NROW,NLAY,KSTP,KPER,TEXT,
     & BUFF,IQ,MXSS,NTSS,NSS,SS,ICBUND,FPRT)
C
C--CHECK IF MAXIMUM NUMBER OF POINT SINKS/SOURCES EXCEEDED.
C--IF SO STOP
  800 WRITE(IOUT,801) NTSS
  801 FORMAT(//1X,'TOTAL NUMBER OF POINT SOURCES/SINKS PRESENT',
     & ' IN THE FLOW MODEL =',I6)
      IF(NTSS.GT.MXSS) THEN
	WRITE(*,802) MXSS
  802	FORMAT(/1X,'ERROR: MAXIMUM NUMBER OF SINKS/SOURCES ALLOWED',
     &	 ' [MXSS] =',I6
     &	 /1X,'INCREASE VALUE OF [MXSS] IN [SSM] PACKAGE INPUT FILE')
C-------EMRL JIG
      call stopfile
C-------EMRL JIG
	STOP
      ENDIF
C
C--IDENTIFY CELLS IN THE VICINITY OF POINT SINKS OR SOURCES
      DO NUM=1,NTSS
	K=SS(1,NUM)
	I=SS(2,NUM)
	J=SS(3,NUM)
	KM1=MAX(K-1,1)
	KP1=MIN(K+1,NLAY)
	DO KK=KM1,KP1
	  IM1=MAX(I-1,1)
	  IP1=MIN(I+1,NROW)
	  DO II=IM1,IP1
	    JM1=MAX(J-1,1)
	    JP1=MIN(J+1,NCOL)
	    DO JJ=JM1,JP1
	      IF(ICBUND(JJ,II,KK,1).EQ.1) ICBUND(JJ,II,KK,1)=1000
	    ENDDO
	  ENDDO
	ENDDO
      ENDDO
C
C--DIVIDE RECH, EVTR, Q_SS BY AQUIFER VOLUME
C--TO GET FLUXES OF SINKS/SOURCES PER UNIT AQUIFER VOLUME.
C--ALSO DETERMINE STEPSIZE WHICH MEETS STABILITY CRITERION
C--FOR SOLVING THE SINK/SOURCE TERM
      DTSSM=1.E30
C
      IF(.NOT.FRCH) GOTO 950
      DO I=1,NROW
	DO J=1,NCOL
	  K=IRCH(J,I)
	  IF(K.EQ.0) CYCLE
	  VOLAQU=DELR(J)*DELC(I)*DH(J,I,K)
	  IF(ICBUND(J,I,K,1).EQ.0.OR.VOLAQU.LE.0) THEN
	    RECH(J,I)=0.
	  ELSE
	    RECH(J,I)=RECH(J,I)/VOLAQU
	  ENDIF
          IF(RECH(J,I).LE.0 .OR. ICBUND(J,I,K,1).EQ.0) CYCLE
	    TM=PRSITY(J,I,K)/RECH(J,I)
	    IF(ABS(TM).LT.DTSSM) THEN
	      DTSSM=ABS(TM)
	      KSSM=K
	      ISSM=I
	      JSSM=J
	    ENDIF
	ENDDO
      ENDDO
C
  950 IF(.NOT.FEVT) GOTO 960
      DO I=1,NROW
	DO J=1,NCOL
	  K=IEVT(J,I)
	  IF(K.EQ.0) CYCLE
	  VOLAQU=DELR(J)*DELC(I)*DH(J,I,K)
	  IF(ICBUND(J,I,K,1).EQ.0.OR.VOLAQU.LE.0) THEN
	    EVTR(J,I)=0
	  ELSE
	    EVTR(J,I)=EVTR(J,I)/VOLAQU
	  ENDIF
          IF(EVTR(J,I).EQ.0 .OR. ICBUND(J,I,K,1).EQ.0) CYCLE
	    TM=PRSITY(J,I,K)/EVTR(J,I)
	    IF(ABS(TM).LT.DTSSM) THEN
	      DTSSM=ABS(TM)
	      KSSM=K
	      ISSM=I
	      JSSM=J
	    ENDIF
	ENDDO
      ENDDO
C
  960 IF(NTSS.LE.0) GOTO 990
      DO NUM=1,NTSS
	K=SS(1,NUM)
	I=SS(2,NUM)
	J=SS(3,NUM)
	VOLAQU=DELR(J)*DELC(I)*DH(J,I,K)
	IF(ICBUND(J,I,K,1).EQ.0.OR.VOLAQU.LE.0) THEN
	  SS(5,NUM)=0
	ELSE
	  SS(5,NUM)=SS(5,NUM)/VOLAQU
	ENDIF
        IF(SS(5,NUM).LE.0 .OR. ICBUND(J,I,K,1).EQ.0) CYCLE
	TM=PRSITY(J,I,K)/SS(5,NUM)
	IF(ABS(TM).LT.DTSSM) THEN
	  DTSSM=ABS(TM)
	  KSSM=K
	  ISSM=I
	  JSSM=J
	ENDIF
      ENDDO
C
C--PRINT INFORMATION ON DTSSM
  990 WRITE(IOUT,1000) DTSSM,KSSM,ISSM,JSSM
 1000 FORMAT(/1X,'MAXIMUM STEPSIZE WHICH MEETS STABILITY CRITERION',
     & ' OF THE SINK & SOURCE TERM'/1X,'=',G11.4,
     & '(WHEN MIN. R.F.=1)  AT K=',I4,', I=',I4,
     & ', J=',I4)
C
C--SYNCHRONIZE ICBUND CONDITIONS OF ALL SPECIES
      IF(NCOMP.EQ.1) GOTO 1999
      DO K=1,NLAY
	DO I=1,NROW
	  DO J=1,NCOL
	    DO INDEX=2,NCOMP
	      IF(ICBUND(J,I,K,INDEX).GE.0) THEN
                ICBUND(J,I,K,INDEX)=IABS(ICBUND(J,I,K,1))
	      ELSEIF(ICBUND(J,I,K,1).EQ.0) THEN
		ICBUND(J,I,K,INDEX)=0
	      ENDIF
	    ENDDO
	  ENDDO
	ENDDO
      ENDDO
C
C--RETURN
 1999 RETURN
      END
C
C
      SUBROUTINE READHQ(INUF,IOUT,NCOL,NROW,NLAY,KSTP,KPER,TEXT,
     & BUFF,FPRT)
C *****************************************************************
C THIS SUBROUTINE READS HEADS AND VOLUMETRIC FLUXES ACROSS CELL
C INTERFACES FROM AN UNFORMATTED FILE SAVED BY THE FLOW MODEL.
C *****************************************************************
C last modified: 06-23-98
C

      IMPLICIT	NONE
      INTEGER	KSTP,KPER,INUF,NCOL,NROW,NLAY,IOUT,IPRTFM,K,I,J,
     &		KKSTP,KKPER,NC,NR,NL
      REAL	BUFF
      CHARACTER TEXT*16,FPRT*1,LABEL*16
      DIMENSION BUFF(NCOL,NROW,NLAY)
	INTEGER*4 RESULT
C
C--WRITE IDENTIFYING INFORMATION
      WRITE(IOUT,1) TEXT,KSTP,KPER,INUF
C
C--READ IDENTIFYING RECORD
      READ(INUF) KKPER,KKSTP,NC,NR,NL,LABEL
C
C--CHECK INTERFACE
      IF(LABEL.NE.TEXT) THEN
	WRITE(*,4) TEXT,LABEL
C-------EMRL JIG
      call stopfile
C-------EMRL JIG
	STOP
      ELSEIF(KKPER.NE.KPER.OR.KKSTP.NE.KSTP) THEN
	WRITE(*,3) KKPER,KKSTP
C-------EMRL JIG
      call stopfile
C-------EMRL JIG
	STOP
      ELSEIF(NC.NE.NCOL.OR.NR.NE.NROW.OR.NL.NE.NLAY) THEN
	WRITE(*,2) NC,NR,NL
C-------EMRL JIG
      call stopfile
C-------EMRL JIG
	STOP
      ENDIF
C
C--READ AN UNFORMATTED RECORD CONTAINING VALUES FOR
C--EACH CELL IN THE GRID
      READ(INUF) (((BUFF(J,I,K),J=1,NCOL),I=1,NROW),K=1,NLAY)
C
C--PRINT OUT INPUT FOR CHECKING IF REQUESTED
      IF(FPRT.NE.'Y'.AND.FPRT.NE.'y') RETURN
      IPRTFM=1
      DO K=1,NLAY
	WRITE(IOUT,50) K
	CALL RPRINT(BUFF(1,1,K),TEXT,
     &	  0,KSTP,KPER,NCOL,NROW,0,IPRTFM,IOUT)
      ENDDO
C
C--PRINT FORMATS
    1 FORMAT(/20X,'"',A16,'" FLOW TERMS FOR TIME STEP',I3,
     & ', STRESS PERIOD',I3,' READ UNFORMATTED ON UNIT',I3
     & /20X,92('-'))
    2 FORMAT(1X,'ERROR: INVALID NUMBER OF COLUMNS, ROWS OR LAYERS',
     & ' IN FLOW-TRANSPORT LINK FILE.'
     & /1X,'NUMBER OF COLUMNS IN FLOW-TRANSPORT LINK FILE =',I5
     & /1X,'NUMBER OF ROWS IN FLOW-TRANSPORT LINK FILE    =',I5,
     & /1X,'NUMBER OF LAYERS FLOW-TRANSPORT LINK FILE     =',I5)
    3 FORMAT(/1X,'ERROR: INVALID NUMBER OF STRESS PERIOD OR TIME STEP',
     &	' IN FLOW-TRANSPORT LINK FILE.'
     & /1X,'NUMBER OF STRESS PERIOD IN FLOW-TRANSPORT LINK FILE =',I3,
     & /1X,'NUMBER OF TIME STEP IN FLOW-TRANSPORT LINK FILE     =',I3)
    4 FORMAT(/1X,'ERROR: INVALID FLOW-TRANSPORT LINK FILE.'/1X,
     & 'NAME OF THE FLOW TERM REQUIRED =',A16/1X,
     & 'NAME OF THE FLOW TERM SAVED IN LINK FILE =',A16)
   10 FORMAT(/44X,'LAYER LOCATION OF ',A16,'FLOW TERM'
     & /44X,43('-'))
   50 FORMAT(/61X,'LAYER ',I3)
C
C--RETURN
      RETURN
      END
C
C
      SUBROUTINE READDS(INUF,IOUT,NCOL,NROW,NLAY,KSTP,KPER,TEXT,
     & BUFF,LOCLAY,FPRT)
C *****************************************************************
C THIS SUBROUTINE READS LOCATIONS AND FLOW RATES OF DIFFUSIVE
C SINK/SOURCE TERMS (RECHARGE AND EVAPOTRANSPIRATION) FROM AN
C UNFORMATTED FILE SAVED BY THE FLOW MODEL.
C *****************************************************************
C last modified: 06-23-98
C

      IMPLICIT	NONE
      INTEGER	KSTP,KPER,INUF,NCOL,NROW,NLAY,IOUT,IPRTFM,I,J,
     &		KKSTP,KKPER,NC,NR,NL,LOCLAY
      REAL	BUFF
      CHARACTER TEXT*16,FPRT*1,LABEL*16
      DIMENSION BUFF(NCOL,NROW),LOCLAY(NCOL,NROW)
	INTEGER*4 RESULT
C
C--WRITE IDENTIFYING INFORMATION
      WRITE(IOUT,1) TEXT,KSTP,KPER,INUF
C
C--READ IDENTIFYING RECORD
      READ(INUF) KKPER,KKSTP,NC,NR,NL,LABEL
C
C--CHECK INTERFACE
      IF(LABEL.NE.TEXT) THEN
	WRITE(*,4) TEXT,LABEL
C-------EMRL JIG
      call stopfile
C-------EMRL JIG
	STOP
      ELSEIF(KKPER.NE.KPER.OR.KKSTP.NE.KSTP) THEN
	WRITE(*,3) KKPER,KKSTP
C-------EMRL JIG
      call stopfile
C-------EMRL JIG
	STOP
      ELSEIF(NC.NE.NCOL.OR.NR.NE.NROW.OR.NL.NE.NLAY) THEN
	WRITE(*,2) NC,NR,NL
C-------EMRL JIG
      call stopfile
C-------EMRL JIG
	STOP
      ENDIF
C
C--READ LAYER LOCATION IF FLOW TERM IS RECHARGE OR E.T.
      READ(INUF) ((LOCLAY(J,I),J=1,NCOL),I=1,NROW)
C
C--READ AN UNFORMATTED RECORD CONTAINING VALUES FOR
C--EACH CELL IN THE GRID
      READ(INUF) ((BUFF(J,I),J=1,NCOL),I=1,NROW)
C
C--PRINT OUT INPUT FOR CHECKING IF REQUESTED
      IF(FPRT.NE.'Y'.AND.FPRT.NE.'y') RETURN
      IPRTFM=1
      CALL RPRINT(BUFF(1,1),TEXT,
     & 0,KSTP,KPER,NCOL,NROW,0,IPRTFM,IOUT)
      IPRTFM=3
      WRITE(IOUT,10)
      CALL IPRINT(LOCLAY(1,1),TEXT,0,KSTP,KPER,NCOL,NROW,
     & 0,IPRTFM,IOUT)
C
C--PRINT FORMATS
    1 FORMAT(/20X,'"',A16,'" FLOW TERMS FOR TIME STEP',I3,
     & ', STRESS PERIOD',I3,' READ UNFORMATTED ON UNIT',I3
     & /20X,92('-'))
    2 FORMAT(1X,'ERROR: INVALID NUMBER OF COLUMNS, ROWS OR LAYERS',
     & ' IN FLOW-TRANSPORT LINK FILE.'
     & /1X,'NUMBER OF COLUMNS IN FLOW-TRANSPORT LINK FILE =',I5
     & /1X,'NUMBER OF ROWS IN FLOW-TRANSPORT LINK FILE    =',I5,
     & /1X,'NUMBER OF LAYERS FLOW-TRANSPORT LINK FILE     =',I5)
    3 FORMAT(/1X,'ERROR: INVALID NUMBER OF STRESS PERIOD OR TIME STEP',
     &	' IN FLOW-TRANSPORT LINK FILE.'
     & /1X,'NUMBER OF STRESS PERIOD IN FLOW-TRANSPORT LINK FILE =',I3,
     & /1X,'NUMBER OF TIME STEP IN FLOW-TRANSPORT LINK FILE     =',I3)
    4 FORMAT(/1X,'ERROR: INVALID FLOW-TRANSPORT LINK FILE.'/1X,
     & 'NAME OF THE FLOW TERM REQUIRED =',A16/1X,
     & 'NAME OF THE FLOW TERM SAVED IN LINK FILE =',A16)
   10 FORMAT(/60X,'LAYER INDEX')
C
C--RETURN
      RETURN
      END
C
C
      SUBROUTINE READPS(INUF,IOUT,NCOL,NROW,NLAY,KSTP,KPER,TEXT,
     & BUFF,IQ,MXSS,NTSS,NSS,SS,ICBUND,FPRT)
C *******************************************************************
C THIS SUBROUTINE READS LOCATIONS AND FLOW RATES OF POINT SINK/SOURCE
C FLOW TERMS FROM AN UNFORMATTED FILE SAVED BY THE FLOW MODEL.
C *******************************************************************
C last modified: 06-23-98
C

      IMPLICIT	NONE
      INTEGER	KSTP,KPER,INUF,NCOL,NROW,NLAY,IOUT,K,I,J,KKSTP,KKPER,
     &		NC,NR,NL,NUM,N,MXSS,NTSS,NSS,ICBUND,IQ,ID
      REAL	BUFF,SS,QSTEMP
      CHARACTER TEXT*16,FPRT*1,LABEL*16
      DIMENSION BUFF(NCOL,NROW,NLAY),ICBUND(NCOL,NROW,NLAY),SS(6,MXSS)
	INTEGER*4 RESULT
C
C--INITIALIZE.
      DO K=1,NLAY
	DO I=1,NROW
	  DO J=1,NCOL
	    BUFF(J,I,K)=0.
	  ENDDO
	ENDDO
      ENDDO
C
C--WRITE IDENTIFYING INFORMATION
      WRITE(IOUT,1) TEXT,KSTP,KPER,INUF
C
C--READ IDENTIFYING RECORD
      READ(INUF) KKPER,KKSTP,NC,NR,NL,LABEL,NUM
C
C--CHECK INTERFACE
      IF(LABEL.NE.TEXT) THEN
	WRITE(*,4) TEXT,LABEL
C-------EMRL JIG
      call stopfile
C-------EMRL JIG
	STOP
      ELSEIF(KKPER.NE.KPER.OR.KKSTP.NE.KSTP) THEN
	WRITE(*,3) KKPER,KKSTP
C-------EMRL JIG
      call stopfile
C-------EMRL JIG
	STOP
      ELSEIF(NC.NE.NCOL.OR.NR.NE.NROW.OR.NL.NE.NLAY) THEN
	WRITE(*,2) NC,NR,NL
C-------EMRL JIG
      call stopfile
C-------EMRL JIG
	STOP
      ENDIF
C
C--RETURN IF NUM=0
      IF(NUM.LE.0) RETURN
C
C--READ AN UNFORMATTED RECORD CONTAINING VALUES FOR
C--EACH POINT SINK OR SOURCE
      DO N=1,NUM
	READ(INUF) K,I,J,QSTEMP
	IF(FPRT.EQ.'Y'.OR.FPRT.EQ.'y') THEN
	  WRITE(IOUT,50) K,I,J,QSTEMP
	ENDIF
	BUFF(J,I,K)=BUFF(J,I,K)+QSTEMP
      ENDDO
C
C--STORE LOCATIONS, FLOW RATES AND TYPES OF
C--THE POINT SINKS OR SOURCES IN ARRAY SS
      DO NUM=1,NSS
	K=SS(1,NUM)
	I=SS(2,NUM)
	J=SS(3,NUM)
	ID=SS(6,NUM)
	IF(ICBUND(J,I,K).EQ.0.OR.BUFF(J,I,K).EQ.0) GOTO 202
	IF(ID.NE.IQ) GOTO 202
	SS(5,NUM)=BUFF(J,I,K)
	IF(ICBUND(J,I,K).LT.0) GOTO 201
	IF(BUFF(J,I,K).LT.0) THEN
	  ICBUND(J,I,K)=1000+IQ
	ELSE
	  ICBUND(J,I,K)=1020+IQ
	ENDIF
  201	BUFF(J,I,K)=0.
  202 ENDDO
      DO K=1,NLAY
	DO I=1,NROW
	  DO J=1,NCOL
	    IF(ICBUND(J,I,K).EQ.0.OR.BUFF(J,I,K).EQ.0) GOTO 203
	    NTSS=NTSS+1
	    IF(NTSS.GT.MXSS) GOTO 203
	    SS(1,NTSS)=K
	    SS(2,NTSS)=I
	    SS(3,NTSS)=J
	    SS(4,NTSS)=0.
	    SS(5,NTSS)=BUFF(J,I,K)
	    SS(6,NTSS)=IQ
	    IF(ICBUND(J,I,K).LT.0) GOTO 203
	    IF(BUFF(J,I,K).LT.0) THEN
	      ICBUND(J,I,K)=1000+IQ
	    ELSE
	      ICBUND(J,I,K)=1020+IQ
	    ENDIF
  203	  ENDDO
	ENDDO
      ENDDO
C
C--PRINT FORMATS
    1 FORMAT(/20X,'"',A16,'" FLOW TERMS FOR TIME STEP',I3,
     & ', STRESS PERIOD',I3,' READ UNFORMATTED ON UNIT',I3
     & /20X,92('-'))
    2 FORMAT(1X,'ERROR: INVALID NUMBER OF COLUMNS, ROWS OR LAYERS',
     & ' IN FLOW-TRANSPORT LINK FILE.'
     & /1X,'NUMBER OF COLUMNS IN FLOW-TRANSPORT LINK FILE =',I5
     & /1X,'NUMBER OF ROWS IN FLOW-TRANSPORT LINK FILE    =',I5,
     & /1X,'NUMBER OF LAYERS FLOW-TRANSPORT LINK FILE     =',I5)
    3 FORMAT(/1X,'ERROR: INVALID NUMBER OF STRESS PERIOD OR TIME STEP',
     &	' IN FLOW-TRANSPORT LINK FILE.'
     & /1X,'NUMBER OF STRESS PERIOD IN FLOW-TRANSPORT LINK FILE =',I3,
     & /1X,'NUMBER OF TIME STEP IN FLOW-TRANSPORT LINK FILE     =',I3)
    4 FORMAT(/1X,'ERROR: INVALID FLOW-TRANSPORT LINK FILE'/1X,
     & 'NAME OF THE FLOW TERM REQUIRED =',A16/1X,
     & 'NAME OF THE FLOW TERM SAVED IN LINK FILE =',A16)
   50 FORMAT(1X,'LAYER',I5,5X,'ROW',I5,5X,'COLUMN',I5,5X,'RATE',G15.7)
C
C--RETURN
      RETURN
      END
C
C
      FUNCTION CREWET(NCOL,NROW,NLAY,CNEW,ICBUND,XBC,YBC,ZBC,
     & JJ,II,KK)
C *****************************************************************
C THIS FUNCTION OBTAINS CONCENTRATION AT A REWET CELL (JJ,II,KK)
C FROM CONCENTRATIONS AT NEIGHBORING NODES WITH INVERSE DISTANCE
C (POWER 2) WEIGHTING .
C *****************************************************************
C last modified: 06-23-98
C
      IMPLICIT	NONE
      INTEGER	NCOL,NROW,NLAY,ICBUND,JJ,II,KK
      REAL	XBC,YBC,ZBC,CTMP,CNEW,CREWET,D2,D2SUM
      DIMENSION ICBUND(NCOL,NROW,NLAY),CNEW(NCOL,NROW,NLAY),
     &		XBC(NCOL),YBC(NROW),ZBC(NCOL,NROW,NLAY)
C
C--INITIALIZE
      D2SUM=0
      CTMP=0
C
C--ACCUMULATE CONCENTRATIONS AT NEIGHBORING NODELS
C--IN THE LAYER DIRECTION
      IF(NLAY.EQ.1) GOTO 10
      IF(KK-1.GT.0) THEN
	IF(ICBUND(JJ,II,KK-1).NE.0) THEN
	  D2=(ZBC(JJ,II,KK)-ZBC(JJ,II,KK-1))**2
	  IF(D2.NE.0) THEN
	    D2SUM=D2SUM+1./D2
	    CTMP=CTMP+CNEW(JJ,II,KK-1)/D2
	  ELSE
	    CTMP=CNEW(JJ,II,KK-1)
	    GOTO 100
	  ENDIF
	ENDIF
      ENDIF
      IF(KK+1.LE.NLAY) THEN
	IF(ICBUND(JJ,II,KK+1).NE.0) THEN
	  D2=(ZBC(JJ,II,KK)-ZBC(JJ,II,KK+1))**2
	  IF(D2.NE.0) THEN
	    D2SUM=D2SUM+1./D2
	    CTMP=CTMP+CNEW(JJ,II,KK+1)/D2
	  ELSE
	    CTMP=CNEW(JJ,II,KK+1)
	    GOTO 100
	  ENDIF
	ENDIF
      ENDIF
C
C--IN THE ROW DIRECTION
   10 IF(NROW.EQ.1) GOTO 20
      IF(II-1.GT.0) THEN
	IF(ICBUND(JJ,II-1,KK).NE.0) THEN
	  D2=(YBC(II)-YBC(II-1))**2
	  IF(D2.NE.0) THEN
	    D2SUM=D2SUM+1./D2
	    CTMP=CTMP+CNEW(JJ,II-1,KK)/D2
	  ELSE
	    CTMP=CNEW(JJ,II-1,KK)
	    GOTO 100
	  ENDIF
	ENDIF
      ENDIF
      IF(II+1.LE.NROW) THEN
	IF(ICBUND(JJ,II+1,KK).NE.0) THEN
	  D2=(YBC(II)-YBC(II+1))**2
	  IF(D2.NE.0) THEN
	    D2SUM=D2SUM+1./D2
	    CTMP=CTMP+CNEW(JJ,II+1,KK)/D2
	  ELSE
	    CTMP=CNEW(JJ,II+1,KK)
	    GOTO 100
	  ENDIF
	ENDIF
      ENDIF
C
C--IN THE COLUMN DIRECTION
   20 IF(NCOL.EQ.1) GOTO 30
      IF(JJ-1.GT.0) THEN
	IF(ICBUND(JJ-1,II,KK).NE.0) THEN
	  D2=(XBC(JJ)-XBC(JJ-1))**2
	  IF(D2.NE.0) THEN
	    D2SUM=D2SUM+1./D2
	    CTMP=CTMP+CNEW(JJ-1,II,KK)/D2
	  ELSE
	    CTMP=CNEW(JJ-1,II,KK)
	    GOTO 100
	  ENDIF
	ENDIF
      ENDIF
      IF(JJ+1.LE.NCOL) THEN
	IF(ICBUND(JJ+1,II,KK).NE.0) THEN
	  D2=(XBC(JJ)-XBC(JJ+1))**2
	  IF(D2.NE.0) THEN
	    D2SUM=D2SUM+1./D2
	    CTMP=CTMP+CNEW(JJ+1,II,KK)/D2
	  ELSE
	    CTMP=CNEW(JJ+1,II,KK)
	    GOTO 100
	  ENDIF
	ENDIF
      ENDIF
C
C--OBTAIN WEIGHTED CONCENTRATION
   30 IF(D2SUM.EQ.0) THEN
	ICBUND(JJ,II,KK)=0
      ELSE
	CTMP=CTMP/D2SUM
      ENDIF
C
C--ASSIGN WEIGHTED CONCENTRATION TO CREWET
  100 CREWET=CTMP
C
C--NORMAL RETURN
      RETURN
      END
